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Nomenclature 


F  inviscid  flux  vector  (axial  direction) 

G  inviscid  flux  vector  (radial  direction) 

K  right  eigenvector  matrix 

L  left  eigenvector  matrix 

R  right-hand-side  residual  vector 

S  source  term  vector 

U  conservative  variable  vector 

A  nozzle  cross  sectional  area 

a  frozen  speed  of  sound 

b  upwind  wave  speed 

E  energy  or  specific  total  energy 

e  specific  energy 

g  degeneracy 

H  specific  total  enthalpy 

h  specific  enthalpy 

i  axial  direction  cell  index 

J  rotational  quantum  number 

j  radial  direction  cell  index 

k  rate  coefficient 

l  cell  side  length 

M  molecular  mass 

m  mass 

p  pressure 

Q  partition  function 

R  specific  gas  constant 

r  radial  coordinate 

T  temperature 

t  time 

u  velocity  component 

v,  w  vibrational  quantum  number 

x  axial  coordinate 

Subscripts 

*  dissociated  state 

oo  free-stream  value 

b  backward 

/  forward 

J  rotational  quantum  number 

L  left  state 

n  cell  side  outward  normal  direction 

R  right  state 

r  radial  direction 

u,  w  vibrational  quantum  number 

x  axial  direction 

Conventions 

c  convection 

D  dissociation 

s  source 

VT  vibrational-translational  energy  exchange 

VV  vibrational- vibrational  energy  exchange 

Symbols 

a  parameter  used  in  the  reconstruction  procedure 
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A  variation 

5  variation 

A  eigenvalue  matrix 

D  cell  volume 

p  density 

a  symmetry  factor 

6  characteristic  rotational  temperature 

Superscripts 

a  atomic  impact  process 

axi  axisymmetric 

chem  chemistry 
m  molecular  impact  process 

v  vibration 

/  formation 

n  time- level 

t  translation 


I.  Introduction 

Possible  applications  of  atmospheric  entry  flows  are  the  design  of  space  vehicle  heat  shields  and  interpre¬ 
tation  of  experimental  data  acquired  in  high  enthalpy  facilities,  such  as  shock  tunnels  and  plasma  torches. 
The  surface  heat  flux  experienced  by  a  spacecraft  during  his  re-entry  represents  the  main  design  parameter. 
Its  value  must  be  accurately  predicted  in  order  to  avoid  mission  failure  and  also  to  reduce  the  margins  to  en¬ 
able  carrying  more  payload.  In  order  to  achieve  this  point,  an  accurate  modeling  of  collisional  and  radiative 
processes  occurring  in  shock  layers  is  needed.  There,  the  flow  experiences  a  sudden  increase  of  temperature, 
pressure  and  density  while  crossing  the  shock  wave.  Atom  and  molecule  collisions  start  occurring  at  higher 
frequencies  and  energies.  This  leads  to  excitation  of  internal  degrees  of  freedom  and  to  molecular  dissocia¬ 
tion  as  well.  If  the  re-entry  speed  is  high  enough,  ionization  may  also  occur.  The  flow  becomes  therefore  a 
partially  ionized  gas,  where  radiative  transitions  take  place.  Their  influence  on  the  medium  dynamics  and 
surface  heat  flux  increases  with  the  increase  of  the  entry  speed.  As  an  example,  for  an  Earth  re-entry  at 
11  km/s,  the  radiative  heat  flux  constitutes  approximately  60%  of  the  total  heating.  The  accurate  description 
of  all  these  phenomena  is  extremely  complex  because  of  nonequilibrium.  This  occurs  in  flow  regions  (such  as 
shock  waves  and  boundary  layers)  where  the  characteristic  time-scales  of  collisional  and  radiative  processes 
become  of  the  same  order  of  magnitude  of  flow  macroscopic  time-scales  (time  needed  for  a  fluid  particle  to 
cross  a  certain  region).  Therefore,  understanding  nonequilibrium  is  crucial  for  a  correct  modeling  of  shock 
layers.  Nonequilibrium  plays  an  important  role  also  in  expanding  flows,  such  as  those  occurring  in  nozzles 
installed  in  high  enthalpy  facilities.  When  heat  flux  and  pressure  measurements  are  performed  on  the  model 
being  tested,  it  is  required  to  know  the  degree  of  nonequilibrium  in  the  free-stream  (generated  by  the  nozzle) 
if  an  extrapolation  to  flight  of  ground  testing  data  is  desired. 

Multi-temperature  models  [T]  were  proposed  as  a  simple  approach  to  overcome  all  difficulties  related  to 
the  modeling  of  nonequilibrium.  They  are  based  on  Boltzmann  distributions  of  the  internal  energy  levels  at 
their  own  temperatures  and  are  valid  only  in  case  of  small  departure  from  equilibrium  [2],  since  the  details 
of  energy  level  dynamics  are  not  taken  into  account.  These  models  are  widely  used  outside  of  their  range  of 
validity  because  of  the  lack  of  more  accurate  physical  models,  but  also  because  of  the  ease  of  implementation 
and  reasonable  amount  of  CPU  time  required  for  3D  simulations. 

State-to-state  or  collisional  models  [SHU],  in  contrast  to  the  multi-temperature  formulation,  treat  each 
energy  level  as  a  separated  pseudo-species.  This  approach  provides  more  flexibility  and  accuracy,  since 
effects  of  non-Boltzmann  distributions  are  accounted  for.  The  price  to  pay  is  that  a  very  large  number  of 
processes  have  to  be  considered  and,  for  each  one,  rate  coefficient  values  have  to  be  obtained  by  means  of 
quantum-mechanical  calculations.  Also,  the  number  of  equations  to  be  solved  rapidly  increases  with  the 
number  of  internal  energy  levels. 

The  main  purpose  of  the  paper  is  to  show  that,  despite  the  higher  computational  time  required,  the 
application  of  collisional  models  to  multi-dimensional  nonequilibrium  flows  has  become  feasible.  The  results 
shown  represent  the  outcome  of  on  going  research  activity  and  are  limited  to  inviscid  flows. 

The  paper  is  structured  as  follows.  Sect.  [TT]  describes  the  physical  model  adopted  (collisional  processes, 
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rate  coefficient  expressions  and  governing  equations  in  conservation  law  form).  The  numerical  method  is 
outlined  in  Sect.  IIIII  Computational  results  are  discussed  in  Sect.  II VI  Conclusions  are  given  in  Sect.  0 

II.  Physical  modeling 

It  is  assumed  that  the  flow  conditions  are  such  that  the  translational  and  rotational  degrees  of  freedom 
of  molecules  are  in  equilibrium.  This  hypothesis  is  not  always  verified  across  shock  waves  [2],  while  it  holds, 
in  general,  for  an  expansion  through  a  converging-diverging  nozzle  in  view  of  the  fast  and  efficient  energy 
exchange  between  rotation  and  translation  at  temperatures  and  pressures  typical  of  nozzles  installed  in  high 
enthalpy  wind  tunnels. 

Chemical  database 

The  starting  point  for  the  development  of  the  vibrational  collisional  (VC)  model  in  use  is  the  chemical 
database  recently  developed  at  NASA  Ames  Research  Center  [3HB],  providing  rate  coefficients  for  the  follow¬ 
ing  processes  involving  the  rovibrational  levels  of  the  nitrogen  molecule  in  its  ground  electronic  state: 

•  Rovibrational  dissociation: 

N2(v,J)  +  Nfe/^*  N  +  N  +  N,  (1) 

•  Rovibrational  excitation: 

t,a 

K  fvJ—tv1  J* 

N2(u,J)+N  N2(t/,J')  +  N,  (2) 

with  v  =  0, . . . ,  vmax,  vf  >  v,  J  =  0, . . . ,  Jmax(v)  and  J'  =  0, . . . ,  Jmax^')-  Rate  coefficients  for  rovibrational 
dissociation  and  excitation  are  computed  by  means  of  the  quasi-classical  trajectory  (QCT)  method  after 
generating  realistic  nuclear  interaction  potentials  by  means  of  quantum-mechanics.  Their  values  are  available 
within  the  range  7500  —  50000  K  (see  J3}jS]  for  more  details). 

The  energy  of  a  rovibrational  level  v  J  is  written  as: 

Evj  =  Ev  4-  A Evj  (3) 

where  Ev  =  Ev o  is  the  vibrational  energy  and  A Evj  =  Evj  —  Ev  is  the  rotational  energy.  The  splitting  in 
eq.  (O  is  not  unique  and  other  choices  are  possible  [lB]. 

Fig.  [T]  shows  the  energy  and  energy  spacing  as  function  of  the  vibrational  quantum  number.  Anhar- 
monicity  effects  can  be  noticed  starting  from  low  quantum  states. 


Figure  1.  Vibrational  energy  levels  and  spacings. 
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Vibrational  collisional  model 

The  VC  model  accounts  for  the  following  processes: 
•  Atomic  impact  dissociation  (Da): 


ka 

N2(t;)  +  N  N  +  N  +  N, 


(4) 


•  Molecular  impact  dissociation  (Dm): 


km 

N2(v)  +  N2  N  +  N  +  N2,  (5) 

k™  v — >-* 


•  Atomic  impact  vibrational-translational  energy  exchange  (VTa): 


ka  / 

N2(v)  +  N  ^  N2(v')+N,  (6) 


•  Molecular  impact  vibrational-translational  energy  exchange  (VTm): 

N2(C+N2fe/^'1N2+N2(t;-l),  (7) 

Kv^v-l 


•  Vibrational- vibrational  energy  exchange  (VV): 

N2(v) +N2(w  —  1)  N2(u  -  1)  +  N2(w),  (8) 


with  v  m  0, . . . ,  umax  (except  for  eq.  ©  where  v  =  1, . . . ,  umax),  vf  «  u  +  1, . . . ,  umax  and  w  =  v,...,  umax. 


Rate  coefficients  and  species  production  rates 

Forward  rate  coefficients  for  Da  and  VTa  processes  in  eqs.  and  @  are  obtained  by  averaging  rovibrational 
rate  coefficients  for  dissociation  and  excitation  and  kJvj^v'jr,  respectively)  over  a  Boltzmann 

distribution  for  rotational  levels  of  each  vibrational  quantum  state  HZ]: 


^  fv— >•* 

^  Tnax(ll) 

=  7T  V  9vJ  exP  ( 

kA  J=0 

v  h’>T 

ka 

"'fv—tv' 

^  ^max(ll) 

=  exp  ( 

^ v  J= 0 

J', 


v  0, . . . ,  umax,  v  <C  v 


(9) 

(10) 


where  k b  is  the  Boltzmann  constant,  gvj  is  the  energy  level  degeneracy  (equal  to  6  ( J  +  1)  for  even  J  and 
3  ( J  +  1)  for  odd  J )  and  Qv  is  the  rotational  partition  function  of  the  vibrational  level  v  [IT] : 


J max  (v) 

Qv=  ^2  9vJ  exp 
j= o 

Backward  rate  coefficients  are  computed  by  means  of  micro-reversibility  m 


VJ  V- 


k\ 


b  v—tv' 


Q^2Qv  f  Ev  —  2F/n  \ 

(ffN  Qn)2  P  \  ^ bT  J  ’ 


4 

Qv 


exp 


Ev  -  Ev>  \ 
fcsT  )  ’ 


e  0)  •  *  *  )  ^max 


(11) 


(12) 

(13) 
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where  E'n  is  the  formation  energy  of  the  nitrogen  atom,  =  12  the  degeneracy  of  its  ground  electronic  state 
(N4S  nuclear  and  electronic  spin  contributions)  and  Q ^  and  Q^2  are  the  translational  partition  functions 
of  the  N  atom  and  N2  molecule,  respectively: 


Qn  ~ 


27rraN&B^  3^2 


Qn2  ~ 


27rraN2&B^  3^2 


(14) 


where  hp  is  the  Planck  constant. 

Rate  coefficients  for  N2  +  N2  inelastic  collisions  are  not  available  in  the  NASA  Ames  database.  It  is  for 
this  reason  that  data  from  literature  are  considered  for  Dm,  VTm  and  VV  processes. 

The  estimation  of  rate  coefficients  for  the  Dm  process  is  performed  by  multi-plying  in  eq.  m  by 

the  Park  factor  /park  =  0.23  [15]: 


kf  v — >•*  —  /Park  kfi 


(15) 


while  the  backward  rate  coefficient  is  obtained  by  means  of  micro-reversibility: 


k™y^*  Qn2Qv  Ev  —  2En  (  N 

=  (16> 

Though  not  rigorous,  this  approach  allows  to  take  into  account  that  N2  is  a  less  effective  collision  partner 
than  N  for  dissociation. 

The  computation  of  forward  rate  coefficients  for  VTm  and  VV  processes  is  performed  by  adapting  the 
existing  set  of  Billing  data  [19]  (based  on  Herzberg  [20]  vibrational  energy  levels)  to  the  vibrational  levels  of 
the  NASA  Ames  database.  The  procedure  is  similar  to  that  in  [2Tj  and  makes  use  of  the  so  called  scaling 
laws  taking  into  account  for  anharmonicity  effects  [22].  Backward  rate  coefficients  are  always  obtained  by 
means  of  micro-reversibility: 


^5  v—tv  —  l 
^ f  v^v  —  1 
L  W-1->W 

^bv^v- 1 


Qv 


exp 


Qv— 1 
Qv  Qw—i 


Ey  Ey—\ 
kpT 


exp  —  - 


Ev  +  Ew   1  —  Ey— 1  —  E, 


kf™zl^w  Qv-iQw~'"  kBT 


V  =  l,...,vn 


W  >  V 


Vfv^v-1 

The  species  production  rates  for  the  N  atom  and  vibrational  levels  of  N2  read: 


cjn  =  2 Mn  Co^a  +  2 Mn  1 


,Dn 


v=G 


v=0 


—  —  Mn2  +  ( 


with  the  partial  contributions  from  all  processes  given  by: 

-'■Da  -  [XN]([Xv]kafv^-[XN]2ktv^*), 
PM]  ([W]  kfv _+*  -  [XN]2  AJV>.)  , 
PM  ([Xv]  kjv^v  -  PM  k%v-+v>)  , 
PM]  ([Xv]  kf -  PM  k^v^v-i 


Urn  a  = 


'Dn 


,VTa 


UJ, 


) 


-VV 


umax 

£  {[Xv]  [Xw-i]  kfv^l-1  +  [Xv-i]  PM  fcMr) ,  «  = 


(17) 

(18) 


(19) 

(20) 

(21) 

(22) 

(23) 

(24) 

(25) 


where  [ X n]  =  and  [Xv\  =  pv/M^2  are  the  concentrations  of  the  N  atom  and  N2  vibrational  levels, 

respectively. 
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Thermodynamics 

The  mixture  density  is  obtained  by  adding  the  contributions  from  the  N  atom  and  N2  vibrational  levels 
p  =  Pn  +  Pn2,  where: 

Pn2  =  ^  Pv  (26) 

v=0 

and  the  static  pressure  follows  from  Dalton’s  law  p  =  p^R^T  p^2R^2T  (where  R n  and  R^2  are  the  specific 
gas  constants).  The  mixture  total  energy  density  is  computed  as: 

3  3  ^max 

pE  =  ~PnRnT  +  —  Pn 2 Rn 2 T  +  ^  pv{&v  +  Aev)  +  pn^n  +  ~^P  (' ux  +  ur)  (27) 

v=0 

where  ux  and  ur  are,  respectively,  the  axial  and  radial  velocity  components,  h £  is  the  specific  formation 
enthalpy  of  the  nitrogen  atom  and  ev  and  Aev  are,  respectively,  the  specific  vibrational  and  rotational  energy 
of  the  vibrational  level  v : 


cv 


Ev 

m^2  ’ 


A  ev(T) 


kBT 2  <9  In Qv 
m^2  dT 


v  0, . . . ,  rmax 


(28) 


The  mixture  total  enthalpy  density  pH  is  obtained  by  adding  the  pressure  to  the  total  energy  density. 

In  order  to  simplify  computations,  the  rigid  rotor  model  is  used  in  eqs.  (fT2|).  (USD,  (USD,  (USD  and  (USD  for 
Qv  and  Aev.  In  this  case  one  has: 


T 

Qv  =  TP  ’  Aev  =  -Rn^T1,  v  =  0, . . . ,  r?max  (29) 

atfY 

where  a  is  the  symmetry  factor  (equal  to  2  for  a  homo-nuclear  diatomic  molecule  such  as  N2),  while  0r  is 
the  characteristic  rotational  temperature  (computed  from  the  rovibrational  energy  levels  as  (Eqi  —  Eqq)  /  kB) . 
The  atomic  nitrogen  degeneracy  g n  in  eq.  m,  must  be  modified  accordingly  and  set  to  4. 

Preliminary  computations  performed  on  quasi  ID  nozzle  and  normal  shock  wave  flows,  have  allowed  to 
verify  that  the  aforementioned  simplification  could  be  applied  for  the  cases  under  investigation. 


Governing  equations 


The  governing  equations  for  2D  axisymmetric  inviscid  nonequilibrium  flows,  written  in  conservation  law 
form,  read: 


dr  U  <9rF  drG  /  . 

^r  +  ^v  +  ^r  =  s  (30) 

where  U,  F  and  G  are,  respectively,  the  conservative  variable  vector  and  the  inviscid  flux  vectors  along  the 
axial  and  radial  directions.  Their  expressions  are  given  by: 


u  = 

[  Pn 

Pv 

P^x  pur 

pE  f, 

(31) 

F  = 

[  Pnux 

Pv^x 

p  +  pul 

piLxUr 

puxH  ]T, 

(32) 

G  = 

[  Pnut 

pvur 

piLxUr 

p  +  pul 

purH  ]T,  v  =  0, . . 

•  i  ^max 

(33) 

The  source  term  S  in  eq.  (I30D  comprises  the  effects  of  the  collisional  processes  previously  described  and  an 
additional  contribution  due  to  the  hypothesis  of  axisymmetric  flow  configuration: 


g  _  gchem  _|_  gaxi 

with  the  chemistry  Schem  and  axisymmetric  Saxl  terms  in  eq.  (Hi  given  by: 

Schem  =  [  rwN  r  U)v  0  0  0  ]T,  V  =  0, . . . ,  vmax, 
Saxi  =  [  0  0  0  p  0  ]T 


(34) 


(35) 

(36) 
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III.  Numerical  method 


Numerical  solutions  to  eq.  (ISO] )  are  obtained  by  means  of  a  2D  multi-block  parallel  Finite  Volume  code 
developed  for  the  application  of  collisional-radiative  models  to  nonequilibrium  flows  [23H25]. 

Governing  equations  are  discretized  in  space  by  means  of  the  Finite  Volume  Method.  This  leads  to  the 
following  semi-discrete  equation  for  the  cell  i,j  (see  Fig.  O  conservative  variable  vector: 


a 


dVt* 


hj 


dt 


D 


(37) 


where  flij  is  the  cell  volrnne  and  R itj  is  the  right-hand-side  residual  made  up  convective  and  source  term 
contributions: 

=  R  ci,j  -I-  Rs ij  (38) 

with: 


Rc  j  7  =  F_  • .  i  •  L,  i  -  4-  F_  •_  i  •  /•_  i  •  +  F_  •  • ,  i  l-  • .  jl  4-  F_  -  -  i  /•  •_  j 

R-sij  =  — 


(39) 

(40) 


In  eq.  (1351.  the  tilde  symbol  (7)  indicates  that  a  numerical  flux  function  is  used  for  the  evaluation  of  convective 
fluxes  at  cell  interfaces  whose  lengths  are  indicated  by  l  (indices  in  eqs.  (1551)  -  (l-TOli  refer  to  Fig.  13  ■ 


i,j  + 1 


Figure  2.  Sketch  of  cell  ordering  for  the  2D  structured  grid  code  in  use. 


Once  the  right-hand-side  residual  split  realized,  the  solution  vector  within  the  cell  i,j  is  updated  at  its 
new  value  at  time  level  n  -h  1  by  adding  contributions  from  the  convection  and  source  term  steps: 


(41) 


In  eq.  (JTT]).  the  two  contributions  to  the  solution  update  are  computed  independently  as  follows.  During 
the  convection  step,  the  presence  of  source  terms  is  neglected  and  the  solution  update  is  computed  explicitly 
according  to: 


Wj  = 


At 


R" 


(42) 


The  convective  time-step  A tcij  in  eq.  ([121)  is  computed  by  means  of  the  usual  CFL-like  condition  [2B]  by 
imposing  for  each  cell  the  maximum  allowable  time-step  allowing  for  local  stability  (local  time-stepping). 

Once  this  done,  the  source  term  step  is  performed  by  neglecting  flux  contributions.  Due  to  the  stiffness 
associated  to  the  source  term,  the  solution  update  is  performed  by  treating  it  implicitly: 


SAX: 


(43) 


where  J  =  dS/d\J  in  eq.  ([11]).  is  the  source  term  Jacobian  matrix.  In  order  to  enhance  stability  and  reduce 
CPU  time,  the  aforementioned  quantity  is  computed  analytically.  The  source  time-step  A tsij  in  eq.  ([IS]) 
should  be  computed  from  stability  considerations  on  the  ODE  system  dXJ/dt  =  S,  and  not  from  the  CFL 
number.  However,  in  the  present  work  the  simplifying  assumption  A =  A tci  j  was  adopted,  since  for 
the  cases  under  investigation  it  did  not  lead  to  a  too  severe  restriction  on  the  time-step. 
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The  approach  here  used  for  time- marching  (known  as  operator  splitting  [26],  has  the  disadvantage,  when 
compared  with  a  fully  implicit  method,  of  having  a  narrower  stability  region.  However,  this  can  be  enlarged 
by  using  multi-stage  time-stepping  schemes  (such  as  the  second  order  two  stage  Adams-Bashforth  method 
[23]  used  in  the  present  work) ,  that  are  known  to  have  a  wider  stability  region  than  the  single-stage  scheme  of 
eq.  dHJ).  Moreover,  the  application  of  the  operator  splitting  approach  does  not  require  the  iterative  solution 
of  large  and  sparse  linear  algebraic  systems,  making  the  solution  update  much  cheaper  than  what  an  implicit 
method  requires. 

The  numerical  flux  in  eq.  03  is  computed  by  means  of  Roe’s  approximate  Riemann  solver  m- 

F„  =  |(F„l  +  Fnii)  -  ^Kra|Ara|L„(Ui{  -  UL)  (44) 

In  eq.  03),  the  normal  physical  flux  Fn  is  computed  according  to  Fn  =  Fn^  +  Gnr  (with  nx  and  nr  being 
the  axial  and  radial  components,  respectively,  of  the  cell  interface  outward  normal  vector),  Kn,  Ln  and  An 
are  the  right  and  left  eigenvector  and  eigenvalue  matrices  projected  along  the  cell  interface  normal  direction, 
while  the  circumflex  accent  (?)  is  used  to  denote  the  Roe-averaged  state  (computed  by  using  the  linearization 
procedure  given  in  [28]).  The  meaning  of  L  and  R  letters  is  obvious  (i.e.  they  correspond  to  the  i  and 
i  +  1  states,  respectively,  if  the  flux  is  being  computed  at  the  interface  i  +  1/2,  j).  Due  to  the  strength  of 
the  shocks  involved,  an  appropriate  entropy  fix  is  required  to  avoid  numerical  instabilities  resulting  from 
even-odd  decoupling.  As  prescribed  by  Quirk  [29].  pressure  disturbances  along  the  shock  front  are  diffused 
away  with  targeted  application  of  the  HLLE  Riemann  solver  [30]: 


F 


n 


b^FnR  —  b  F  nL 

b+  —  b~ 


b+b~ 
b+  —  b~ 


(U«  -  U L) 


(45) 


where  the  upwind  wave  speeds  b  are  given  by: 

(46) 

(47) 


b+  =  max(0,  VnR  +  aR,  Vn  +  a), 
b~  =  min(0,  VnL  -  ah,Vn  -  a) 


In  eqs.  03  -  03,  Vn  is  the  normal  velocity  computed  as  Vn  =  uxnx  +  urnr ,  while  a  represents  the  frozen 
speed  of  sound  of  the  mixture. 

High-order  spatial  resolution  is  achieved  by  means  of  parabolic  interpolation  of  the  left  and  right  states 
at  the  cell  interfaces  via: 

U  l  —  -(2U^_i  j  +  5LJ —  Ui+1J)  (48) 

The  use  of  eq.  ([48])  makes  the  scheme  3rd-order  accurate  in  space  m •  Due  to  the  limiting  required  for 
strong  nonlinear  waves,  the  reconstructed  left  and  right  states  must  be  modified  according  to: 


UL  <-  median(UL,  U itj,  UMP)  (49) 

with: 

UMP  =  U ij  +  mimnodtUi+ij  -  Ui?i,  a(Vij  -  Ui-ij)}  (50) 

being  the  monotonicity-preserving  (MP)  limit  [32]  and  minmod(a,  b)  =  ^(sgn(a)+sgn(5))min(abs(a),  abs(b)). 
For  the  present  work,  the  parameter  a  has  been  set  to  2.  The  right  state  can  be  easily  found  from  symmetry. 
Due  to  the  non-linearity  of  the  limiter,  the  reconstruction  is  performed  on  the  characteristic  variables, 
requiring  a  full  diagonalization  of  the  governing  equations.  More  details  are  given  in  [23]. 


IV.  Computational  results 

The  VC  model  described  in  Sec.  [II] has  been  applied  for  the  computation  of  the  following  2D  axisymmetric 
inviscid  flows: 

•  Supersonic  flow  through  a  converging-diverging  nozzle, 

•  Flow  over  a  sphere. 
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For  the  first  test-case,  the  nozzle  of  the  NASA  Ames  EAST  facility  has  been  chosen  (see  Fig.  [3J),  while  for 
the  second  one  a  sphere  of  radius  R  =  0.5  m  has  been  considered.  In  order  to  get  more  physical  insight 
from  both  simulations,  a  vibrational  temperature  T*  has  been  extracted  from  the  vibrational  energy  level 
populations  by  solving  the  following  algebraic  non-linear  equation: 

%ax  %ax 

Y.  nvEv  Ev  exp 

u=0 _  _  v=0 _ 

vmax  Umax  / 

Y  £  exp  ( 

v=0  v=0  \ 


at  each  mesh  cell. 


Figure  3.  Normalized  area  distribution  of  the  EAST  facility  nozzle. 


Nozzle  flow 

The  flow  is  supposed  to  be  in  equilibrium  at  the  nozzle  inlet.  The  reservoir  total  pressure  and  temperature 
are,  respectively,  po  =  100  atm  and  To  =  5600  K  and  correspond  to  actual  operative  conditions  of  the  EAST 
facility.  For  the  latters,  experimental  data  on  vibrational  energy  level  populations  (measured  by  means  of 
Raman  spectroscopy)  are  available  [33]. 

Computations  were  run  by  using  the  4-block  mesh  shown  in  Fig.  [5]  Each  block  consists  of  a  50x50 
cell  grid.  In  order  to  assess  the  solution  accuracy  on  that  grid,  a  convergence  study  has  been  performed  by 
doubling  the  number  of  cells  (for  each  block)  in  both  radial  and  axial  directions.  The  solutions  obtained  by 
using  the  original  and  finer  grids  did  not  show  appreciable  departure  from  each  other,  hence  justifying  the 
use  of  the  original  grid. 

Figs.  [5] -ED  show  the  velocity  magnitude  and  temperature  (translational  and  vibrational)  fields  together 
with  the  distribution  of  the  same  quantities  along  the  nozzle  axis.  The  flow  remains  in  thermal  equilibrium 
for  all  the  converging  part  of  the  nozzle,  as  can  be  appreciated  from  the  temperature  distribution.  Once 
it  crosses  the  throat,  the  expansion  becomes  significant  and  nonequilibrium  effects  start  to  appear.  This  is 
caused  by  the  reduction  of  the  flow  macroscopic  time-scales  due  to  the  velocity  increase.  The  translational  and 
vibrational  temperatures  start  deviating  from  each  other,  with  the  latter  becoming  frozen  around  cm. 

Downstream  this  location,  the  flow  further  expands  as  if  the  inelastic  collisional  processes  taken  into  account 
(see  Sec.  |TT|  were  not  occurring. 

Figs.  [7]  show  the  mass  fractions  of  the  N  atom,  N2  molecule  and  the  first  9  vibrational  levels  of  the  latter. 
The  degree  of  recombination  is  quite  low.  This  is  due  to  the  fact  that  the  flow  in  the  reservoir  is  weakly 
dissociated  as  an  effect  of  the  high  and  low  values  of  total  pressure  and  temperature,  respectively.  Despite 
that,  an  accurate  prediction  of  the  aforementioned  quantity  is  always  needed,  since  a  too  approximated 
prediction  of  the  flow  outlet  conditions  can  lead  to  an  erroneous  interpretation  of  experimental  data  acquired 
in  high-enthalpy  wind  tunnels.  From  the  same  figures,  one  can  see  that  the  relative  population  of  the  ground 
state  increases,  while  those  of  excited  states  experience  depletion. 
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0.02- 


Figure  4.  Multi-block  mesh  used  for  the  EAST  facility  nozzle  (4  blocks  —  200x50  cells). 


This  does  not  mean  that  the  recombination  occurs  mainly  in  the  ground  state.  As  a  matter  of  fact, 
the  observed  behavior  of  the  excited  vibrational  levels  is  a  consequence  of  the  concurrent  action  of  Da,  Dm 
processes,  leading  to  the  formation  of  vibrationally  excited  molecules,  and  VTa,  VTm  and  W  inelastic 
processes  that,  as  an  effect  of  the  cooling  the  flow  is  undergoing,  lead  to  an  increase  of  the  population  of 
low-lying  levels. 


0.05  -| 


(a)  Velocity. 


(b)  Velocity  along  nozzle  axis. 


Figure  5.  Velocity  field. 


The  usage  of  a  VC  model  allows  for  a  detailed  investigation  of  the  energy  level  dynamics.  Moreover,  since 
no  hypothesis  is  done  on  the  population  distribution,  it  is  possible  to  check  the  approximation  introduced 
by  classical  multi-temperature  models  (where  a  Boltzmann  distribution  at  its  own  temperature  is  supposed 
to  exist  for  each  species  internal  degree  of  freedom).  Fig.  HI  shows  the  evolution  of  the  distribution  along  the 
nozzle  axis.  At  the  inlet,  where  equilibrium  is  imposed,  all  levels  lie  along  a  straight  line.  When  the  flow 
begins  to  expand  in  the  converging  portion  of  the  nozzle,  the  shape  of  the  distribution  remains  the  same.  Only 
the  slope  changes  as  an  effect  of  the  cooling.  A  small  overpopulation  starts  to  appear  (for  high-lying  levels)  at 
the  throat  and  this  becomes  more  and  more  pronounced  when  the  flow  becomes  supersonic.  At  the  outlet,  the 
ground  state  and  low-lying  levels  are  aligned  along  a  straight  line  (whose  slope  is  proportional  to  the  inverse 
of  the  vibrational  temperature  T*  previously  shown)  while  high-lying  levels  experience  overpopulation.  The 
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latter  is  an  indication  that  lumping  all  the  levels  within  a  Boltzmann  distribution  (as  done  within  the  context 
of  multi-temperature  models)  is  not  possible  (for  the  reservoir  conditions  adopted  here). 


0.02 


1  I  1 

0.02 


1 — I — r- 

0.04 


0.06 


0.08 


(a)  Temperatures  (upper  part  T  -  lower  part  Tv  ). 

Figure  6.  Translational  and  vibrational  temperature  fields. 


Figure  7.  Mass  fraction  evolution  along  the  nozzle  axis. 


The  computed  population  distributions  have  been  compared  against  experimental  data.  These  have  been 
obtained  by  means  of  Raman  spectroscopy  [33]  at  the  locations  provided  in  Table  U 


Table  1.  Locations  of  experimental  data  acquisition. 


1  2  3 

x  [cm]  -0.6  2.4  5.4 


Fig.  EH  compares  the  computed  and  experimental  normalized  population  distributions  for  the  first  10  vi¬ 
brational  levels  of  the  N2  molecule.  The  agreement  is  excellent  for  the  first  location  (converging  portion 
of  the  nozzle),  where  the  flow  is  close  to  equilibrium.  For  the  second  and  third  location,  the  agreement  is 
good  for  the  ground  state  and  the  first  four  excited  levels,  while  discrepancies  start  to  appear  for  the  higher 
states.  In  this  zone  (where  measurement  errors  are  more  significant),  the  computational  model  adopted 
predicts  a  Boltzmann  distribution  at  the  temperature  Tv,  while  an  overpopulation  appears  in  experimental 
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data.  The  latter  could  be  due  to  the  presence  of  some  nonequilibrium  at  the  nozzle  inlet.  This  possible 
source  of  disagreement  is  totally  neglected  in  the  computations,  since  equilibrium  is  assumed  at  nozzle  inlet. 
Moreover,  the  actual  flow  in  the  facility  is  characterized  by  the  presence  of  viscous  effects  and  unsteadiness. 
None  of  these  factors  are  taken  into  account  within  the  computational  model  adopted.  In  view  of  that,  the 
comparison  between  computational  and  experimental  data  can  be  considered  satisfactory. 


Figure  8.  Population  evolution  along  the  nozzle  axis  ( - x  —  —2.5  cm, - x  =  —0.6  cm, - x  =  0  cm, - x  —  2.4  cm, 

- - - x  —  5.4cm,  •  •  •  •  x  —  8.3cm). 


Figure  9.  Comparison  between  computed  and  experimental  normalized  vibrational  energy  level  population  distribu¬ 
tions. 


13  of  E| 


American  Institute  of  Aeronautics  and  Astronautics 


Flow  over  a  sphere 

For  the  computation  of  the  inviscid  2D  axisymmetric  flow  over  a  sphere,  the  free- stream  conditions  provided 
in  Table  [2]  have  been  adopted: 


Table  2.  Free-stream  conditions. 

Poo  [Pa]  Tqq  [K]  Fqq  [m/s] 

13.3  300  9,000 


A  12-block  mesh  (see  Fig.  [TO),  with  each  block  consisting  of  a  90x30  cell  grid,  has  been  used. 


Figure  10.  Multi-block  mesh  used  for  the  sphere  (12  blocks  —  90  x  360  cells) . 


Figs.  HU-EU  show  the  temperature  distribution  and  their  evolution  along  stagnation  line,  respectively. 
At  the  shock  location,  the  vibrational  temperature  is  frozen  and  maintains  its  free-stream  value.  Behind  the 
shock,  excitation  occurs  at  the  beginning  mainly  through  VTm  processes. 


0.8 


0.6 


0.4 
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X 

(a)  Translational  temperature. 


X 

(b)  Vibrational  temperature. 


Figure  11.  Translational  and  vibrational  temperature  fields. 
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Figure  12.  Temperature  evolution  along  the  stagnation  line. 


Once  molecules  are  sufficiently  excited,  they  start  to  dissociate.  In  this  phase,  as  soon  N  is  formed  because 
of  dissociation,  VTa  processes  start  to  occur  as  well.  The  collisional  excitation  leads  to  an  initial  increase 
of  the  vibrational  temperature  since  medium  and  high-lying  are  being  populated.  On  the  other  hand,  the 
dissociation  causes  a  depletion  of  the  population  of  vibrational  levels  and  tends  to  decrease  the  vibrational 
temperature  value.  After  an  initial  stage  where  molecules  are  being  excited,  the  dissociation  becomes  the 
dominant  mechanism  and  causes  the  vibrational  temperature  (as  shown  in  Fig.  [T2])  to  reach  a  maximum 
and  then  relax  and  assume  the  same  value  as  the  translational  temperature.  Since  the  flow  is  inviscid, 
no  boundary  layer  is  present,  and  the  wall  temper ature  assumes  the  value  corresponding  to  equilibrium 
conditions  downstream  the  shock.  It  is  interesting  to  notice  that,  as  opposed  to  the  results  obtained  by  means 
of  multi-temperature  models,  the  vibrational  temperature  does  not  become  higher  than  the  translational  one. 
The  aforementioned  behavior  is  not  physical  and  is  due  to  the  use  of  too  approximated  formulations  for  the 
chemistry- vibration  coupling  source  term  pQ.  When  one  uses  a  VC  model  (like  in  the  present  case),  there  is 
no  need  for  this  macroscopic  source  term,  since  energy  levels  are  treated  as  pseudo-species  and  the  chemistry- 
vibration  coupling  is  consistently  and  intrinsically  contained  in  the  rate  coefficients  for  dissociation  from  all 
vibrational  levels. 

The  dynamics  of  the  macroscopic  dissociation  of  the  N2  molecule  can  also  be  appreciated  from  Fig.  [T51 
showing  the  evolution  along  the  stagnation  line  of  the  mass  fractions  of  the  N  atom,  N2  molecule  and  the 
first  10  vibrational  levels  of  the  latter.  These  provide  a  confirmation  to  what  was  mentioned  before,  that  is, 
macroscopic  dissociation  proceeds  through  an  initial  stage  of  excitation  followed  by  a  second  phase  where 
dissociation  from  all  vibrational  levels  becomes  the  dominant  mechanism. 


Figure  13.  Mass  fraction  evolution  along  the  stagnation  line. 
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Fig.  [TU  shows  the  evolution  of  the  vibrational  energy  level  population  between  the  shock  and  wall 
locations.  The  initial  shape  corresponds  to  a  Boltzmann  distribution  at  the  free-stream  temperature.  As 
soon  as  the  flow  crosses  the  shock,  collisional  excitation  (through  single  and  multi-quantum  transitions)  starts 
populating  medium  and  high-lying  levels.  At  this  stage  the  population  of  vibrational  levels  strongly  deviates 
from  a  Boltzmann  distribution.  After  that,  N2  molecules  start  dissociating  and  the  population  of  vibrational 
levels  approaches  the  Boltzmann  distribution  corresponding  to  the  post-shock  equilibrium  state.  The  present 
analysis,  despite  the  simplicity  of  the  flow  field  being  computed,  shows  that  macroscopic  dissociation  in  a 
shock  layer  flow  does  not  occur  through  a  sequence  of  Boltzmann  distributions  of  vibrational  levels. 


Figure  14.  Population  evolution  along  the  stagnation  line  ( -  x  =  —0.55  m, - x  —  —0.53  m,  —  •  —  x  =  —0.525  m, 

—  .  .  —  x  —  —0.52m, - - - x  —  —0.515cm,  ....  x  =  0m). 


V.  Conclusions 

A  vibrational  collisional  model  has  been  developed  and  applied  for  investigating  nonequilibrium  effects 
in  2D  axisymmetric  inviscid  flows.  Applications  have  considered  the  supersonic  flow  through  a  converging- 
diverging  nozzle  (EAST  facility  of  NASA  Ames  Research  Center)  and  the  flow  past  a  sphere.  In  both 
cases,  the  detailed  analysis  of  the  vibrational  energy  level  dynamics  has  shown  that,  both  macroscopic  re¬ 
combination  and  dissociation  do  not  proceed  through  a  sequence  of  Boltzmann  distributions  at  their  own 
temperatures  (as  supposed  within  the  context  of  multi- temperature  models). 

The  computed  vibrational  energy  level  population  distributions  for  the  EAST  facility  nozzle  have  been 
compared  against  experimental  data.  Despite  the  simplifying  assumption  introduced  in  the  modeling  phase 
(inviscid  and  steady  flow)  a  fair  agreement  has  been  observed. 

Future  work  will  focus  on  the  implementation  of  viscous  effects  and  comparison  with  results  obtained  by 
means  of  multi-temperature  models. 
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